

%% load estimation results
load('parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
lambdaStar = theta((size(thetaStar,1)-5):end,use);
UStar = payoffs(thetaStar,X,Z,YY,reb,mpow);
clear K l0 srngSV SV theta theta0 time use Y exitn Z X K;


%% load everything of interest
load('sim_eq_final.mat')
[DMES,intervEq,EqPayoffs] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);

%% preallocate
B = size(DMES,1);
PRpeace = NaN(B, N, 4);
Enumenter = NaN(B, N, 4);
Enumentercond = NaN(B, N, 4);
IntervenersData = NaN(B,N,5);
IntervenersBR = NaN(B,N,5);
IntervenersMP = NaN(B,N,5);
IntervenersMW = NaN(B,N,5);

for n=1:N
    for b=1:NS
        if ~isempty(intervEq{b,n})
             
            EQ=intervEq{b,n};
            spne=intervEq{b,n}(1,:)>=V_R(1,b,n);
            EQ(1,spne)=0;
            EQ(:,~spne)=repmat([1;zeros(size(EQ,1)-1,1)],1,sum(~spne));
            ExpFighters = EQ'*sum(YY,2);
            ExpFighterscon = intervEq{b,n}(2:end,:)'*sum(YY(2:end,:),2);
              
            Ref = sortrows([ExpFighters ExpFighterscon (1:size(EQ,2))']);
            
            sgpayoffs = intervEq{b,n}(2:end,:)'*[V_R(2:end,b,n) V_US(2:end,b,n) V_UK(2:end,b,n)...
                                               V_France(2:end,b,n) V_Russia(2:end,b,n) V_China(2:end,b,n)];
            
            % cf_1: simulated data
            PrA = profprob(1:size(YY,1), DMES{b,n}, intervEq{b,n}, V_R(1,b,n), lambdaStar);
            PRpeace(b,n,1) = PrA(1,1);
            Enumenter(b,n,1) = PrA' * sum(YY(:,2:end),2);
            peq = exp(DMES{b,n} * lambdaStar);
            peq = peq / sum(peq);
            Enumentercond(b,n,1) = sum(sum((peq' .* intervEq{b,n}(2:end,:)) .* sum(YY(2:end,2:end),2)));
            IntervenersData(b,n,:) = sum((intervEq{b,n}(2:end,:) * peq) .* YY(2:end,2:end));
            
            % cf_2: max peace
            eq = EQ(:, Ref(1,3));
            neqsg = intervEq{b,n}(:, Ref(1,3));
            PRpeace(b,n,2) = unique(eq(1,:));
            Enumenter(b,n,2) = unique(eq' * sum(YY(:,2:end),2));
            Enumentercond(b,n,2) =  unique(neqsg(2:end,:)' * sum(YY(2:end,2:end),2));
            IntervenersMP(b,n,:) = sum(neqsg(2:end,:) .* YY(2:end,2:end));
        
            % cf_3: max war
            eq = EQ(:, Ref(end,3));
            neqsg = intervEq{b,n}(:, Ref(end,3));
            PRpeace(b,n,3) = unique(eq(1,:));
            Enumenter(b,n,3) = unique(eq' * sum(YY(:,2:end),2));
            Enumentercond(b,n,3) =  unique(neqsg(2:end,:)' * sum(YY(2:end,2:end),2));
            IntervenersMW(b,n,:) = sum(neqsg(2:end,:) .* YY(2:end,2:end));
            
            % cf_4: best for the rebels
            RefR = sortrows([EqPayoffs{b,n}(:,1) sgpayoffs(:,1) (1:size(EQ,2))']);
            eq = EQ(:, RefR(end,3));
            neqsg = intervEq{b,n}(:, RefR(end,3));
            PRpeace(b,n,4) = unique(eq(1,:));
            Enumenter(b,n,4) = unique(eq' * sum(YY(:,2:end),2));
            Enumentercond(b,n,4) =  unique(neqsg(2:end,:)' * sum(YY(2:end,2:end),2));
            IntervenersBR(b,n,:) = sum(neqsg(2:end,:) .* YY(2:end,2:end));    
        end
    end
end


%%

clearvars -except max_peace max_war best_R best_US best_UK best_FR best_RU... 
    best_CH PRpeace Enumenter Enumentercond...
    IntervenersBR IntervenersData IntervenersMP IntervenersMW N

load('data.mat')
data.Properties.VariableNames{1,1} = 'ccode';
Enumentercond(isinf(Enumentercond))  = NaN;
meanPRpeace=NaN(N,size(PRpeace,3));
meanEnumenter = NaN(N,size(Enumenter,3));
meanEnumentercond = NaN(N,size(Enumentercond,3));


for i = 1:4
   meanPRpeace(:,i) = nanmean(PRpeace(:,:,i));
   meanEnumenter(:,i) = nanmean(Enumenter(:,:,i));
   meanEnumentercond(:,i) = nanmean(Enumentercond(:,:,i));
end


cf_names = {'data', 'maxPeace', 'maxWar',  'bestRebel'};
meanPRpeace = array2table(meanPRpeace, 'VariableNames', cf_names);
meanPRpeace = [data(:,1:2), meanPRpeace];

meanEnumenter = array2table(meanEnumenter, 'VariableNames', cf_names);
meanEnumenter = [data(:,1:2), meanEnumenter];

meanEnumentercond = array2table(meanEnumentercond, 'VariableNames', cf_names);
meanEnumentercond = [data(:,1:2), meanEnumentercond];

writetable(meanPRpeace, 'selection_prpeace_replication.csv')

